home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 November / macformat-018.iso / Utility Spectacular / Utilities / Calc / xmath.h < prev    next >
Encoding:
C/C++ Source or Header  |  1993-10-10  |  14.7 KB  |  440 lines  |  [TEXT/????]
  1. /*
  2.  * Copyright (c) 1992 David I. Bell
  3.  * Permission is granted to use, distribute, or modify this source,
  4.  * provided that this copyright notice remains intact.
  5.  *
  6.  * Data structure declarations for extended precision arithmetic.
  7.  * The assumption made is that a long is 32 bits and shorts are 16 bits,
  8.  * and longs must be addressible on word boundaries.
  9.  */
  10.  
  11. #include "alloc.h"
  12.  
  13. #include "have_stdlib.h"
  14. #ifdef HAVE_STDLIB_H
  15. #include <stdlib.h>
  16. #endif
  17.  
  18.  
  19. #ifndef    NULL
  20. #define    NULL    0
  21. #endif
  22.  
  23. /*#define ALLOCTEST 1*/
  24.  
  25. #ifndef ALLOCTEST
  26. # if defined(UNIX_MALLOC)
  27. #  define freeh(p) { if ((p != _zeroval_) && (p != _oneval_)) free((void *)p); }
  28. # else
  29. #  define freeh(p) ((p == _zeroval_) || (p == _oneval_) || free(p))
  30. # endif
  31. #endif
  32.  
  33. typedef    short FLAG;            /* small value (e.g. comparison) */
  34. typedef unsigned short BOOL;        /* TRUE or FALSE value */
  35.  
  36. #if !defined(TRUE)
  37. #define    TRUE    ((BOOL) 1)            /* booleans */
  38. #endif
  39. #if !defined(FALSE)
  40. #define    FALSE    ((BOOL) 0)
  41. #endif
  42.  
  43.  
  44. /*
  45.  * NOTE: FULL must be twice the storage size of a HALF
  46.  *     LEN storage size must be <= FULL storage size
  47.  */
  48. typedef unsigned short HALF;        /* unit of number storage */
  49. typedef unsigned long FULL;        /* double unit of number storage */
  50. typedef long LEN;            /* unit of length storage */
  51.  
  52. #define BASE    ((FULL) 65536)        /* base for calculations (2^16) */
  53. #define BASE1    ((FULL) (BASE - 1))    /* one less than base */
  54. #define BASEB    16            /* number of bits in base */
  55. #define    BASEDIG    5            /* number of digits in base */
  56. #define    MAXHALF    ((FULL) 0x7fff)        /* largest positive half value */
  57. #define    MAXFULL    ((FULL) 0x7fffffff)    /* largest positive full value */
  58. #define    TOPHALF    ((FULL) 0x8000)        /* highest bit in half value */
  59. #define    TOPFULL    ((FULL) 0x80000000)    /* highest bit in full value */
  60. #define MAXLEN    ((LEN)    0x7fffffff)    /* longest value allowed */
  61. #define    MAXREDC    5            /* number of entries in REDC cache */
  62. #define    SQ_ALG2    20            /* size for alternative squaring */
  63. #define    MUL_ALG2 20            /* size for alternative multiply */
  64. #define    POW_ALG2 40            /* size for using REDC for powers */
  65. #define    REDC_ALG2 50            /* size for using alternative REDC */
  66.  
  67.  
  68.  
  69. typedef union {
  70.     FULL    ivalue;
  71.     struct {
  72.         HALF Svalue1;
  73.         HALF Svalue2;
  74.     } sis;
  75. } SIUNION;
  76.  
  77.  
  78. #ifdef LITTLE_ENDIAN
  79.  
  80. #define silow    sis.Svalue1    /* low order half of full value */
  81. #define sihigh    sis.Svalue2    /* high order half of full value */
  82.  
  83. #else
  84.  
  85. #define silow    sis.Svalue2    /* low order half of full value */
  86. #define sihigh    sis.Svalue1    /* high order half of full value */
  87.  
  88. #endif
  89.  
  90.  
  91. typedef struct {
  92.     HALF    *v;        /* pointer to array of values */
  93.     LEN    len;        /* number of values in array */
  94.     BOOL    sign;        /* sign, nonzero is negative */
  95. } ZVALUE;
  96.  
  97.  
  98.  
  99. /*
  100.  * Function prototypes for integer math routines.
  101.  */
  102. #if defined(__STDC__)
  103. #define proto(a) a
  104. #else
  105. #define proto(a) ()
  106. #endif
  107.  
  108. extern HALF * alloc proto((LEN len));
  109. #ifdef    ALLOCTEST
  110. extern void freeh proto((HALF *));
  111. #endif
  112.  
  113.  
  114. extern long iigcd proto((long i1, long i2));
  115. extern void itoz proto((long i, ZVALUE * res));
  116. extern long ztoi proto((ZVALUE z));
  117. extern void zcopy proto((ZVALUE z, ZVALUE * res));
  118. extern void zadd proto((ZVALUE z1, ZVALUE z2, ZVALUE * res));
  119. extern void zsub proto((ZVALUE z1, ZVALUE z2, ZVALUE * res));
  120. extern void zmul proto((ZVALUE z1, ZVALUE z2, ZVALUE * res));
  121. extern void zsquare proto((ZVALUE z, ZVALUE * res));
  122. extern void zreduce proto((ZVALUE z1, ZVALUE z2,
  123.     ZVALUE * z1res, ZVALUE * z2res));
  124. extern void zdiv proto((ZVALUE z1, ZVALUE z2,
  125.     ZVALUE * res, ZVALUE * rem));
  126. extern void zquo proto((ZVALUE z1, ZVALUE z2, ZVALUE * res));
  127. extern void zmod proto((ZVALUE z1, ZVALUE z2, ZVALUE * rem));
  128. extern BOOL zdivides proto((ZVALUE z1, ZVALUE z2));
  129. extern void zor proto((ZVALUE z1, ZVALUE z2, ZVALUE * res));
  130. extern void zand proto((ZVALUE z1, ZVALUE z2, ZVALUE * res));
  131. extern void zxor proto((ZVALUE z1, ZVALUE z2, ZVALUE * res));
  132. extern void zshift proto((ZVALUE z, long n, ZVALUE * res));
  133. extern long zlowbit proto((ZVALUE z));
  134. extern long zhighbit proto((ZVALUE z));
  135. extern BOOL zisset proto((ZVALUE z, long n));
  136. extern BOOL zisonebit proto((ZVALUE z));
  137. extern BOOL zisallbits proto((ZVALUE z));
  138. extern void zbitvalue proto((long n, ZVALUE * res));
  139. extern FLAG ztest proto((ZVALUE z));
  140. extern FLAG zrel proto((ZVALUE z1, ZVALUE z2));
  141. extern BOOL zcmp proto((ZVALUE z1, ZVALUE z2));
  142. extern void trim proto((ZVALUE * z));
  143. extern void shiftr proto((ZVALUE z, long n));
  144. extern void shiftl proto((ZVALUE z, long n));
  145. extern void zfact proto((ZVALUE z, ZVALUE * dest));
  146. extern void zpfact proto((ZVALUE z, ZVALUE * dest));
  147. extern void zlcmfact proto((ZVALUE z, ZVALUE * dest));
  148. extern void zperm proto((ZVALUE z1, ZVALUE z2, ZVALUE * res));
  149. extern void zcomb proto((ZVALUE z1, ZVALUE z2, ZVALUE * res));
  150. extern BOOL zprimetest proto((ZVALUE z, long count));
  151. extern FLAG zjacobi proto((ZVALUE z1, ZVALUE z2));
  152. extern void zfib proto((ZVALUE z, ZVALUE * res));
  153. extern void zpowi proto((ZVALUE z1, ZVALUE z2, ZVALUE * res));
  154. extern void ztenpow proto((long power, ZVALUE * res));
  155. extern void zpowermod proto((ZVALUE z1, ZVALUE z2,
  156.     ZVALUE z3, ZVALUE * res));
  157. extern BOOL zmodinv proto((ZVALUE z1, ZVALUE z2, ZVALUE * res));
  158. extern void zgcd proto((ZVALUE z1, ZVALUE z2, ZVALUE * res));
  159. extern void zlcm proto((ZVALUE z1, ZVALUE z2, ZVALUE * res));
  160. extern BOOL zrelprime proto((ZVALUE z1, ZVALUE z2));
  161. extern long zlog proto((ZVALUE z1, ZVALUE z2));
  162. extern long zlog10 proto((ZVALUE z));
  163. extern long zdivcount proto((ZVALUE z1, ZVALUE z2));
  164. extern long zfacrem proto((ZVALUE z1, ZVALUE z2, ZVALUE * rem));
  165. extern void zgcdrem proto((ZVALUE z1, ZVALUE z2, ZVALUE * res));
  166. extern long zlowfactor proto((ZVALUE z, long count));
  167. extern long zdigits proto((ZVALUE z1));
  168. extern FLAG zdigit proto((ZVALUE z1, long n));
  169. extern BOOL zsqrt proto((ZVALUE z1, ZVALUE * dest));
  170. extern void zroot proto((ZVALUE z1, ZVALUE z2, ZVALUE * dest));
  171. extern BOOL zissquare proto((ZVALUE z));
  172. extern void zmuli proto((ZVALUE z, long n, ZVALUE *res));
  173. extern long zmodi proto((ZVALUE z, long n));
  174. extern long zdivi proto((ZVALUE z, long n, ZVALUE * res));
  175. extern HALF *zalloctemp proto((LEN len));
  176.  
  177. #if 0
  178. extern void zapprox proto((ZVALUE z1, ZVALUE z2, ZVALUE* res1, ZVALUE* res2));
  179. #endif
  180.  
  181.  
  182. /*
  183.  * Modulo arithmetic definitions.
  184.  * Structure holding state of REDC initialization.
  185.  * Multiple instances of this structure can be used allowing
  186.  * calculations with more than one modulus at the same time.
  187.  * Len of zero means the structure is not initialized.
  188.  */
  189. typedef    struct {
  190.     LEN len;        /* number of words in binary modulus */
  191.     ZVALUE mod;        /* modulus REDC is computing with */
  192.     ZVALUE inv;        /* inverse of modulus in binary modulus */
  193.     ZVALUE one;        /* REDC format for the number 1 */
  194. } REDC;
  195.  
  196. #if 0
  197. extern void zmulmod proto((ZVALUE z1, ZVALUE z2, ZVALUE z3, ZVALUE *res));
  198. extern void zsquaremod proto((ZVALUE z1, ZVALUE z2, ZVALUE *res));
  199. extern void zsubmod proto((ZVALUE z1, ZVALUE z2, ZVALUE z3, ZVALUE *res));
  200. #endif
  201. extern void zminmod proto((ZVALUE z1, ZVALUE z2, ZVALUE *res));
  202. extern BOOL zcmpmod proto((ZVALUE z1, ZVALUE z2, ZVALUE z3));
  203. extern REDC *zredcalloc proto((ZVALUE z1));
  204. extern void zredcfree proto((REDC *rp));
  205. extern void zredcencode proto((REDC *rp, ZVALUE z1, ZVALUE *res));
  206. extern void zredcdecode proto((REDC *rp, ZVALUE z1, ZVALUE *res));
  207. extern void zredcmul proto((REDC *rp, ZVALUE z1, ZVALUE z2, ZVALUE *res));
  208. extern void zredcsquare proto((REDC *rp, ZVALUE z1, ZVALUE *res));
  209. extern void zredcpower proto((REDC *rp, ZVALUE z1, ZVALUE z2, ZVALUE *res));
  210.  
  211.  
  212. /*
  213.  * Rational arithmetic definitions.
  214.  */
  215. typedef struct {
  216.     ZVALUE num, den;
  217.     long links;
  218. } NUMBER;
  219.  
  220. extern NUMBER *qadd(), *qsub(), *qmul(), *qdiv(), *qquo(), *qmod(), *qcomb();
  221. extern NUMBER *qsquare(), *qgcd(), *qlcm(), *qmin(), *qmax(), *qand(), *qor();
  222. extern NUMBER *qxor(), *qpowermod(), *qpowi(), *qpower(), *qneg(), *qsign();
  223. extern NUMBER *qfact(), *qpfact(), *qsqrt(), *qshift(), *qminv();
  224. extern NUMBER *qint(), *qfrac(), *qnum(), *qden(), *qinv(), *qabs(), *qroot();
  225. extern NUMBER *qfacrem(), *qcopy(), *atoq(), *itoq(), *iitoq();
  226. extern NUMBER *qperm(), *qgcdrem(), *qtrunc(), *qround(), *qalloc();
  227. extern NUMBER *qlowfactor(), *qfib(), *qcfappr(), *qcos(), *qsin(), *qexp();
  228. extern NUMBER *qscale(), *qln(), *qbtrunc(), *qbround(), *qisqrt();
  229. extern NUMBER *qtan(), *qacos(), *qasin(), *qatan(), *qatan2(), *qjacobi();
  230. extern NUMBER *qinc(), *qdec(), *qhypot(), *qcosh(), *qsinh(), *qtanh();
  231. extern NUMBER *qacosh(), *qasinh(), *qatanh(), *qlegtoleg(), *qiroot();
  232. extern NUMBER *qpi(), *qbappr(), *qdivi(), *qlcmfact(), *qminmod();
  233. extern NUMBER *qredcin(), *qredcout(), *qredcmul(), *qredcsquare();
  234. extern NUMBER *qredcpower();
  235. extern BOOL qcmp(), qcmpi(), qprimetest(), qissquare();
  236. extern BOOL qisset(), qcmpmod(), qquomod();
  237. extern FLAG qrel(), qreli(), qnear(), qdigit();
  238. extern long qtoi(), qprecision(), qplaces(), qdigits();
  239. extern long qilog2(), qilog10(), qparse();
  240. extern void qfreenum();
  241. extern void qprintnum();
  242. extern void setepsilon();
  243.  
  244. #if 0
  245. extern NUMBER *qbitvalue(), *qmuli(), *qmulmod(), *qsquaremod();
  246. extern NUMBER *qaddmod(), *qsubmod(), *qreadval(), *qnegmod();
  247. extern BOOL qbittest();
  248. extern FLAG qtest();
  249. #endif
  250.  
  251. #ifdef CODE
  252. extern NUMBER *qaddi();
  253. #endif
  254.  
  255.  
  256. /*
  257.  * Complex arithmetic definitions.
  258.  */
  259. typedef struct {
  260.     NUMBER *real;        /* real part of number */
  261.     NUMBER *imag;        /* imaginary part of number */
  262.     long links;        /* link count */
  263. } COMPLEX;
  264.  
  265. extern COMPLEX *cadd(), *csub(), *cmul(), *cdiv(), *csquare();
  266. extern COMPLEX *cneg(), *cinv();
  267. extern COMPLEX *comalloc(), *caddq(), *csubq(), *cmulq(), *cdivq();
  268. extern COMPLEX *cpowi(), *csqrt(), *cscale(), *cshift(), *cround();
  269. extern COMPLEX *cbround(), *cint(), *cfrac(), *croot(), *cexp(), *cln();
  270. extern COMPLEX *ccos(), *csin(), *cpolar(), *cpower(), *cmodq(), *cquoq();
  271. extern void comfree(), comprint();
  272. extern BOOL ccmp();
  273. extern void cprintfr();
  274.  
  275. #if 0
  276. extern COMPLEX *cconj(), *creal(), *cimag(), *qqtoc();
  277. #endif
  278.  
  279.  
  280. /*
  281.  * macro expansions to speed this thing up
  282.  */
  283. #define iseven(z)    (!(*(z).v & 01))
  284. #define isodd(z)    (*(z).v & 01)
  285. #define iszero(z)    ((*(z).v == 0) && ((z).len == 1))
  286. #define isneg(z)    ((z).sign)
  287. #define ispos(z)    (((z).sign == 0) && (*(z).v || ((z).len > 1)))
  288. #define isunit(z)    ((*(z).v == 1) && ((z).len == 1))
  289. #define isone(z)    ((*(z).v == 1) && ((z).len == 1) && !(z).sign)
  290. #define isnegone(z)    ((*(z).v == 1) && ((z).len == 1) && (z).sign)
  291. #define istwo(z)    ((*(z).v == 2) && ((z).len == 1) && !(z).sign)
  292. #define isleone(z)    ((*(z).v <= 1) && ((z).len == 1))
  293. #define istiny(z)    ((z).len == 1)
  294. #define issmall(z)    (((z).len < 2) || (((z).len == 2) && (((short)(z).v[1]) >= 0)))
  295. #define isbig(z)    (((z).len > 2) || (((z).len == 2) && (((short)(z).v[1]) < 0)))
  296. #define z1tol(z)    ((long)((z).v[0]))
  297. #define z2tol(z)    (((long)((z).v[0])) + \
  298.                 (((long)((z).v[1] & MAXHALF)) << BASEB))
  299.  
  300. #define qiszero(q)    (iszero((q)->num))
  301. #define qisneg(q)    (isneg((q)->num))
  302. #define qispos(q)    (ispos((q)->num))
  303. #define qisint(q)    (isunit((q)->den))
  304. #define qisfrac(q)    (!isunit((q)->den))
  305. #define qisunit(q)    (isunit((q)->num) && isunit((q)->den))
  306. #define qisone(q)    (isone((q)->num) && isunit((q)->den))
  307. #define qisnegone(q)    (isnegone((q)->num) && isunit((q)->den))
  308. #define qistwo(q)    (istwo((q)->num) && isunit((q)->den))
  309. #define qiseven(q)    (isunit((q)->den) && iseven((q)->num))
  310. #define qisodd(q)    (isunit((q)->den) && isodd((q)->num))
  311. #define qistwopower(q)    (isunit((q)->den) && zistwopower((q)->num))
  312. #define qhighbit(q)    (zhighbit((q)->num))
  313. #define qlowbit(q)    (zlowbit((q)->num))
  314. #define qdivides(q1, q2)    (zdivides((q1)->num, (q2)->num))
  315. #define qdivcount(q1, q2)    (zdivcount((q1)->num, (q2)->num))
  316. #define qilog(q1, q2)    (zlog((q1)->num, (q2)->num))
  317. #define qlink(q)    ((q)->links++, (q))
  318.  
  319. #define qfree(q)    {if (--((q)->links) <= 0) qfreenum(q);}
  320. #define quicktrim(z)    {if (((z).len > 1) && ((z).v[(z).len-1] == 0)) (z).len--;}
  321.  
  322. #define cisreal(c)    (qiszero((c)->imag))
  323. #define cisimag(c)    (qiszero((c)->real) && !cisreal(c))
  324. #define ciszero(c)    (cisreal(c) && qiszero((c)->real))
  325. #define cisone(c)    (cisreal(c) && qisone((c)->real))
  326. #define cisnegone(c)    (cisreal(c) && qisnegone((c)->real))
  327. #define cisrunit(c)    (cisreal(c) && qisunit((c)->real))
  328. #define cisiunit(c)    (qiszero((c)->real) && qisunit((c)->imag))
  329. #define cistwo(c)    (cisreal(c) && qistwo((c)->real))
  330. #define cisint(c)    (qisint((c)->real) && qisint((c)->imag))
  331. #define ciseven(c)    (qiseven((c)->real) && qiseven((c)->imag))
  332. #define cisodd(c)    (qisodd((c)->real) || qisodd((c)->imag))
  333. #define clink(c)    ((c)->links++, (c))
  334.  
  335. #define    clearval(z)    memset((z).v, 0, (z).len * sizeof(HALF))
  336. #define    copyval(z1, z2)    memcpy((z2).v, (z1).v, (z1).len * sizeof(HALF))
  337.  
  338.  
  339. /*
  340.  * Flags for qparse calls
  341.  */
  342. #define QPF_SLASH    0x1    /* allow slash for fractional number */
  343. #define QPF_IMAG    0x2    /* allow trailing 'i' for imaginary number */
  344.  
  345.  
  346. /*
  347.  * Output modes for numeric displays.
  348.  */
  349. #define MODE_DEFAULT    0
  350. #define MODE_FRAC    1
  351. #define MODE_INT    2
  352. #define MODE_REAL    3
  353. #define MODE_EXP    4
  354. #define MODE_HEX    5
  355. #define MODE_OCTAL    6
  356. #define MODE_BINARY    7
  357. #define MODE_MAX    7
  358.  
  359. #define MODE_INITIAL    MODE_REAL
  360.  
  361.  
  362. /*
  363.  * Output routines for either FILE handles or strings.
  364.  */
  365. extern void math_chr(), math_str(), math_flush();
  366. extern void divertio(), cleardiversions(), setfp();
  367. extern char *getdivertedio();
  368. extern void setmode();        /* set output mode for numeric output */
  369. extern void setdigits();    /* set # of digits for float or exp output */
  370.  
  371. #ifdef VARARGS
  372. extern void math_fmt();
  373. #else
  374. # ifdef __STDC__
  375. extern void math_fmt(char *, ...);
  376. # else
  377. extern void math_fmt();
  378. # endif
  379. #endif
  380.  
  381. /*
  382.  * Print a formatted string containing arbitrary numbers, similar to printf.
  383.  */
  384. #ifdef VARARGS
  385. extern void qprintf();
  386. #else
  387. # ifdef __STDC__
  388. extern void qprintf(char *, ...);
  389. # else
  390. extern void qprintf();
  391. # endif
  392. #endif
  393.  
  394. /*
  395.  * constants used often by the arithmetic routines
  396.  */
  397. extern HALF _zeroval_[], _oneval_[], _twoval_[], _tenval_[];
  398. extern ZVALUE _zero_, _one_, _ten_;
  399. extern NUMBER _qzero_, _qone_, _qnegone_, _qonehalf_;
  400. extern COMPLEX _czero_, _cone_;
  401.  
  402. #if 0
  403. extern NUMBER _conei_;
  404. #endif
  405.  
  406. extern BOOL _sinisneg_;        /* whether sin(x) < 0 (set by cos(x)) */
  407. extern BOOL _math_abort_;    /* nonzero to abort calculations */
  408. extern long _epsilonprec_;    /* binary precision of epsilon */
  409. extern NUMBER *_epsilon_;    /* default error for real functions */
  410. extern ZVALUE _tenpowers_[32];    /* table of 10^2^n */
  411. extern long _outdigits_;    /* current output digits for float or exp */
  412. extern int _outmode_;        /* current output mode */
  413. extern LEN _mul2_;        /* size of number to use multiply algorithm 2 */
  414. extern LEN _sq2_;        /* size of number to use square algorithm 2 */
  415. extern LEN _pow2_;        /* size of modulus to use REDC for powers */
  416. extern LEN _redc2_;        /* size of modulus to use REDC algorithm 2 */
  417. extern HALF *bitmask;        /* bit rotation, norm 0 */
  418.  
  419. #if 0
  420. extern char *_mallocptr_;    /* pointer for malloc calls */
  421. #endif
  422.  
  423. /*
  424.  * misc function declarations - most to keep lint happy
  425.  */
  426. extern void initmasks();    /* init the bitmask rotation arrays */
  427.  
  428. #ifdef VARARGS
  429. void error();
  430. #else
  431. # ifdef __STDC__
  432. void error(char *, ...);
  433. # else
  434. void error();
  435. # endif
  436. #endif
  437.  
  438.  
  439. /* END CODE */
  440.